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Abstract: Strongly-coupled fermionic systems can support a variety of low-energy phe- 
nomena, giving rise to collective condensation, symmetry breaking and a rich phase struc- 
ture. We explore the potential of worldline Monte Carlo methods for analyzing the effec- 
tive action of fermionic systems at large flavor number A'f, using the Gross-Neveu model 
as an example. Since the worldline Monte Carlo approach does not require a discretized 
spacetime, fermion doubling problems are absent, and chiral symmetry can manifestly be 
maintained. As a particular advantage, fluctuations in general inhomogeneous conden- 
sates can conveniently be dealt with analytically or numerically, while the renormalization 
can always be uniquely performed analytically. We also critically examine the limitations 
of a straightforward implementation of the algorithms, identifying potential convergence 
problems in the presence of fermionic zero modes as well as in the high-density region. 
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1. Introduction 



The elaboration of the QCD phase diagram as a function of temperature and baryon density 
represents a great challenge in modern particle and nuclear physics. It has been studied 
experimentally at the RHIC collider at Brookhaven and at CERN. It has also attracted 
intense theoretical investigation using computer simulations in the context of lattice gauge 
theory and random matrix theory . A significant computational obstacle is the so-called 
sign problem, for which a number of approximate techniques have been developed: (i) a 
Taylor expansion with respect to the baryon chemical potential |Q, |3|, (ii) imaginary 
values for the chemical potential /x fs], |^ ; (iii) overlap enhancing techniques . Despite 
these impressive successes, our knowledge of the QCD phase diagram is still rather limited. 
Other approaches include the study of QCD-like theories such as 2-colour QCD |lO|, 
11], and high density quark matter in a colour superconducting state [12, |l^, |l^, |l5|, 
16]. Of particular interest is the identification of inhomogeneous phases, analogous to the 
LOFF states of condensed matter systems [|^, |l^, 19, 20]. This idea has been studied 
deeply for Skyrme models ]|2l|, 23, 24], and for four-fermion models ]|2^, 26, | 



. Lower 

dimensional models also provide insight into chiral dynamics, and early lattice studies of 
lower-dimensional models such as the Gross Neveu model in the limit of many flavours have 
shown that inhomogeneous background fields are key to understand the dense phase ]27]. 
Recently, this model has attracted a lot of interest since it was found analytically that the 
high-density state of fermion matter forms an inhomogeneous "baryon crystal" |2^, 29, 30, 
31]. This analytic crystalline phase has been confirmed by a lattice analysis and such 



studies have led to several new insights concerning inhomogeneous phases [|33| , p4| ]. 

Good chiral properties of the fermion action is of central importance for an investigation 
of quark matter at intermediate densities, since the high-density transition is driven by 
chiral dynamics. Unfortunately, lattice fermion actions necessarily suffer from the fermion 
doubling problem as firstly pointed out by Nielsen and Ninomiya ]35]. Nowadays, staggered 
fermions p^ , domain wall fermions ]37] or Neuberger fermions |3^, which are an explicit 
realization of the Ginsparg- Wilson relation |3£], are widely used in numerical simulations. 
Despite these advanced formulations and great numerical efforts, it turns out still to be 
cumbersome to achieve good chiral properties, such as a sufficiently small pion mass. 

Here we propose to use the worldline approach to study interacting fermion systems. 
Since the worldline approach to the quark determinant does not use a lattice discretization 
of space-time, it circumvents many of these significant difficulties. Here, we will argue 
that the prospects of the worldline approach are (i) exact chiral symmetry but yet a fully 
numerical approach, (ii) analytic renormalization and (iii) a clear description of Fermi 
surface effects. The worldline method is a string-inspired approach to quantum field theory; 
see [40, ^ for reviews. It was further developed into a viable tool for an efficient calculation 
of functional determinants for arbitrary background fields [42[. Subsequently, worldline 



numerics has enjoyed a wide span of applications ranging from the Casimir effect [43, 44[ 



and fermion induced quantum interactions [45[ to the description of pair production in 



inhomogeneous fields [46, 47[. A worldline lattice formulation has also been presented in 

Hi. 
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In this paper we apply a worldline Monte Carlo approach to interacting fermion models, 
concentrating on the (1 + l)-dimensional Gross-Neveu model 5C]. We review these 



models in Section 2, and describe the worldline representation of the effective action. The 
treatment of inhomogeneous condensates is presented in Section 3, and the extension to 
finite temperature and chemical potential is explained in Section 3, along with results of the 
computations. Section 5 contains our conclusions, and several appendices contain relevant 
technical details. 

2. Worldline representation of the effective action 
2.1 Continuum formulation 

In this work, we mainly concentrate on the the D = 1 + 1 dimensional Gross-Neveu (GN) 
model defined by the Euclidean fermionic action |49, pH] 



Sf = I + (2.1) 

where A'f denotes the number of flavors. The GN model has a discrete chiral symmetry 
under ip — > 75'i/'- Our conventions for the Dirac algebra in D = 1 -|- 1 are 71 = ia^, 72 = io"i, 
75 = o'2i {7^;7i^} = —"^S^iu, but most of the equations below generalize straightforwardly 
to higher dimensions. With the aid of a Hubbard-Stratonovich transformation, we can 
rewrite this model in terms of a mixed fermionic-bosonic action 



2g2 



(2.2) 



with the bosonic scalar (j{x) which transforms as a — > —a under the discrete chiral symme- 
try. In this work, we actually confine ourselves strictly to the large-A^f limit of this model 
which is equivalent to a semiclassical approximation of the a sector. Integrating out the 
remaining fermion fluctuations, results in the large- A^f effective action 

T[a]= [ d^x^a'^ - NfTi\n{-a-ia), ^ 00. (2.3) 
J 

The resulting Dirac operator V = ii^ — a has a 75 hermiticity, = 752575, implying that 
the fermionic contribution to F is real. In this case, we can use the identity TrlniP = 
iTrlnPPt to obtain 

r[a] = [ d^x^a^ - ^Trln(-92 + - i^a), N{ ^ 00, (2.4) 
J 25^ 2 

now involving a fluctuation operator of Laplace type. 

A variety of fermionic models of a similar type exist. An example for a model with a 
continuous chiral symmetry is given by the Thirring model which is discussed in appendix 

Here, we concentrate on the Gross-Neveu model. The contribution from the fermion 
fluctuations to the effective action can be written in propertime form, 

Ar[a] := -^Trln(-52 +cj2 -i^CT) (2.5) 

^ r ^^^^-i~^^+^^-^^<^)T^ (2.6) 



2 -/i/A2 T 
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where we have used a propertime UV cutoff A for the sake of definiteness; any other 
regularization scheme can similarly be implemented. Interpreting the fluctuation operator 
in the exponent as a quantum mechanical Hamiltonian H, we can introduce the path 
integral representation of the propertime evolution operator exp(— iJT), 



^rl-l = f(I^'^/r^ / f--'^"^**'"»--'-^*"'*W. (2.7) 

1/A2 x(0)=x{T) 

containing the spin factor 



1 



T 



$[0-] = — tr^ T^exp ( i / dr {^a) ] , {2.i 



and a are functions of x(r). The normalization of the path integral in Eq. (2.7) has been 
chosen such that f Vx exp(— f x'^/A) = 1. We also introduced the dimensionality of the 
Dirac algebra dy, e.g., dy = 2^/^ in even dimensions, such that $[0] = 1. 

2.2 Propertime discretization 

For a given background of arbitrarily varying a and V fields, the analytic computation of the 
effective action T is formidable. For a numerical computation, the worldline representation 
is highly useful, since the computation of the path integral in Eqs. ( |2.7D or ( A.5| ) can 
be done efficiently with Monte Carlo methods. Most importantly, the algorithm can be 
formulated for any background a in full generality. 

An estimate for the path integral over an operator 0[x] is given by an average over a 
finite ensemble of closed worldlines ("loops") X£{t), where £ = 1, . . . ,nL. 

/T 1 
VxO[x]e'yo dT±Hr) = ^ (2.9) 

Since the path integral is over all spacetime, it can be decomposed into a spacetime integral 
over a "center of mass" (CM) and a path integral over worldlines fixed to this common 
center of mass, 

{0[x]) = J d^xcM(OW).eM (2-10) 

This latter ensemble of mass-centered worldlines xi is generated according to a Gauf5ian 
velocity distribution 

Pcm[x] = 5 (^xcM - £dTxir)^ e-llo'i^^'(r), (2.11) 

where the S function implements the center-of-mass constraint. The center-of-mass integra- 
tion relates the effective action with a corresponding action density, i.e., the Lagrangian, 
r = J d^xcwC-ixcu), but this decomposition of the worldlines is not unique. Alterna- 
tively, the worldlines in the ensemble can be fixed to have one common point (CP), say at 
r = 0, 

Pcp[x] = 6 (xcp - x{t = 0)) e~yo'i^^'i^)^ (2.12) 
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such that (©[x]) = f d^xcp{0[x])xcp- The associated Lagrangian C,T = f d^xcpC{xcp), 
is not identical to the one related to the center-of-mass definition, but agrees with it up to 
a total derivative term . Below, we use both prescriptions: whereas the CM worldlines 
have some computational advantages due to the constraint J x{t) = 0, analytical results 
are often related to the representation involving CP worldlines. 

The finite ensemble of worldline loops still has an infinite set of degrees of freedom. We 
reduce these to finitely many by a discretization of the propertime coordinate r, resulting 
in an approximation of each worldline by a finite set of points per loop, Xi := x(rj), Tj = 
Ti/N, i = 0, . . . , N, where the points Xi=o = Xi=N are identified. Replacing the propertime 
derivative in the worldline action by a nearest-neighbor difference, i;(rj) — > N{xi — Xi-i)/T, 
leads to a GauBian action for the Xi which can be diagonalized including the 6 constraint. 
This allows for an ab-initio generation of a random worldline distribution; for concrete 
algorithms, see [43, 51]. A discretization in terms of Fourier modes is also possible and 
will be used in Sect. S.C; it shows a slightly slower convergence towards the propertime 
continuum limit ||4^ . 

Whereas the finiteness of the number of worldlines introduces an error which is measur- 
able by the statistical error of the expectation values, the propertime discretization leads 
to a systematic error which has to be studied by approaching the continuum limit, i.e., by 
increasing the value of A^. 

It should be stressed that only the auxiliary propertime is discretized, whereas the 
spacetime coordinates remain continuous, Xi G (at least to machine accuracy). There- 
fore, there is neither a fermion doubling problem nor is the algorithm restricted to certain 
values of the flavor number A'f ; in fact, A'^f can even be a continuous number. 

Another advantage of the worldline numerical approach is given by the fact that chiral 
symmetries can manifestly be preserved. For instance in the case of the Thirring model 
[discussed in more detail in Appendix^, the SU(iVf)RX SU(A'f)L remains unaffected also 
in the discretized version of the worldline approach. This is because the derivation of 
the worldline representation is completely done in the continuum, each step preserving 
the chiral properties of the Dirac operator. The final worldline expression ( [A.5| ) hence is 
manifestly chirally invariant. Discretizing the propertime does not spoil this invariance; 
moreover, spacetime remains continuous even after propertime discretization. Even if we 
decided to discretize the worldline integrals by constraining them to form chains of links 
on a spacetime lattice, chiral symmetry would not be violated, since such a discretization 
would not correspond to a lattice Dirac operator. We conclude that the numerical worldline 
Monte Carlo estimate for the effective action respects chiral invariance, has no doubling 
problem and can be used for any number of flavors A'f . 



2.3 Renormalization 

Renormalization of these fermionic models in the large- A^f limit can be done in the standard 
fashion, using, e.g., Coleman- Weinberg renormalization conditions to fix physical param- 
eters. The large-A^f limit even facilitates renormalization in any dimension D, requiring, 
of course, an increasing number of counter terms corresponding to physical parameters for 
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increasing D} 

In practice, however, renormalization can become an issue if the computation is done 
numerically, since divergent pieces have to be separated from the physical quantities in a 
well-controlled manner. Therefore, we briefly outline the renormalization process in the 
following in order to demonstrate that the renormalization can be done analytically first 
in our formalism, such that numerics has to deal only with the remaining finite parts. For 
simplicity, we confine ourselves to D = 1 + 1 in the Gross-Neveu model but the same 
strategy can be applied also in the case of other fermionic models. 

The unrenormalized large- A/'f effective action in the Gross-Neveu model reads 

oo 

^ 1/A2 

where we have already included the free-field subtraction, such that r[(T = 0] = 0, and 
labeled the bare coupling with a subscript A. A logarithmic cutoff dependence is generated 
from the small-T limit of the propertime integrand. The small-T expansion of the worldline 
expectation value, which corresponds to the heat-kernel expansion, yields 

{e-fod^'''q>[a] - 1)^ = -a\x)T + 0{T^), (2.14) 

where only the term linear in T supports a log A divergence. This divergence corresponds 
to the local counter term proportional to the bare action and can be absorbed by a cou- 
pling renormalization. We fix the physical coupling at the renormalization scale // by a 
renormalization condition of Coleman- Weinberg type: 



6a^ 



(2.15) 



Inserting Eq. ( 2.13| ), this relates the bare to the renormalized coupling. 



1 1 / a'^ 



2 +^ {C + 2 + \n^] , (2.16) 

where C denotes the Euler-Mascheroni constant. Also the /3 function can be read off, Pg2 = 
fidi^tg{fi)'^ = —^g{fj,)^. The finite parts are specific to our propertime cutoff regularization; 
of course, any other desired regularization method can also be used, such that the final 
numerical results can be mapped to any given regularization scheme. 

The renormalization scale and the coupling can be traded for an RG invariant mass 
scale, manifesting dimensional transmutation. A convenient choice is given by the minimum 
m of the effective action which is directly related to the induced fermion mass, 

6r 

5a 



= =^ m = /ue^ . (2.17) 

=m,A— +00 



^Beyond the large- A^f limit, the models are perturbatively renormalizable in D = 1 + 1. For instance in 
D = 2 + 1, they are even nonperturbatively renormalizable owing to the occurrence of a non-Gaufiian fixed 
point, as can be proven, e.g., to all order in a 1/Nf expansion |53, psj]. 
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Expressing the bare coupling in terms of the renormahzed couphng and successively in 
terms of the minimum m, leads us to the final renormalized form of the effective action, 



T[a] 



47r 



oo 



a' 







r2 



•x 9I 



.-rr?T _ 1) _ 1^ 



(2.18) 



which holds for any given background a{x). Note that the worldline expectation value 
that is to be computed numerically is finite. Together with the subtraction terms, also 
the propertime integral is finite in the limit A ^ 00, which we have performed already in 
Eq. ( 2.18| ). For constant o", we have (e~-^o '^'^'^^ ^\a\ix = e~"^'^ , and the effective action can 
be calculated analytically, 



r[a = const.] = I Sx 
47r 



a 

In^ -1 



(2.19) 



exhibiting a nontrivial minimum at cr = m giving rise to chiral symmetry breaking and 
fermion mass generation. 

For the present work, the representation Eq. ( |2.18| ) is not only conceptually satisfactory, 
but also useful from a practical viewpoint, since the minimum m which sets the physical 
scale is under control analytically, cf. Eq. ( 2.19| ). For the remainder of this subsection, we 
briefly consider an alternative but equally valid choice of the scale. 

The disadvantage of the above scale fixing lies in the fact that m is determined from the 
minimum of the full effective action, which happens to be at a =constant in the large- A'f 
limit. In the general case, e.g., beyond large A'^f or in a fully nonperturbative calculation, 
the true minimum of F may not be known (or at least relatable to the renormalized cou- 
pling) beforehand. Still, a suitable subtraction scale is required for the computation. An 
alternative choice is, for instance, given by the average field ^'^ = j (fxa{x)'^, where 0,2 
is the 1 + 1 dimensional spacetime volume. With this choice, the large-A'f renormalized 
effective action becomes 



T\a] 



47r 



d^x 



a 



In 



CJ 



1 + 



dT 







(2.20) 



where m is still defined via Eq. (2.17). This choice of the subtraction scale has been used 
in |^3| . Beyond large A'^f or fully nonperturbatively, the first log term will be replaced by a 
more complicated effective action, but the fermion loop contribution can still be computed 
in this manner and remains finite. 



3. Spatially inhomogeneous a condensates 
3.1 Worldline formalism 

A number of nontrivial tests of the preceding general worldline approach can already be 
performed within the Gross-Neveu model in Z) = 1 + 1 which we will exclusively consider 
in the following. For the case of static but spatially varying a condensates, a number of 
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exact results can be derived, and these serve as stringent tests and controls of our general 
worldline approach. Special attention will be paid to bound states, especially zero modes 
and near-zero modes, as these pose the greatest numerical challenges. 

For a static but spatially varying a{x), the renormalized action in the 1 + 1 dimensional 
Gross-Neveu model at large iVf boils down to 

(3.1) 

where the prime denotes the derivative with respect to the spatial coordinate, a' = da/dx. 
It is convenient to combine exponential and hyperbolic function in the expectation value 
by defining 

V±{x) :=a^{x)±a'{x) , (3.2) 



and 



^ = yJ d''x{jr+ + C.) , (3.3) 



where 

:^ -V. . [ % ((.-/-.)^ - - 1) - l) . (3.4, 

So far, we have not specified the details of the worldline ensemble. For a comparison with 
analytical results, common-point (CP) loops are particularly useful, since then we can 
write the associated expectation value as a quantum mechanical transition amplitude in 
imaginary time, 

-^'^^^^) = ^I'^rf 1^)1 ^1^^ = ^{x\e-^^-\x) ^ V4^K^{x;T) . (3.5) 
' occF [x\e ^ ^'\x) 

Here we have introduced the Hamiltonian H± = —df + V±, and the heat kernel K±{x;T), 
with X = xcp being the common point. We have normalized the heat kernel such that 
it agrees with the standard conventions of ID quantum mechanical transition amplitudes. 
For example, with a homogeneous condensate a{x) = m, we have V± = m^, and 

The transition amplitude, or heat kernel, in turn can be written as a sum over the eigen 
modes ipn of the Hamiltonian H±: 

K±{x-T) ^ (x|e-^^±|x) = |V'n(^)|'e-^^S, (3.7) 

n 

where the eigenvalues of the Hamiltonian are denoted by E"^. This spectral decomposition 
fixes the large-T behavior of the worldline expectation value. If the Hamiltonian H± has 
no zero mode, then all terms of the sum are exponentially damped and the expectation 
value vanishes for large T. By contrast, if a zero mode exists and the spectrum is discrete 
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at its lower end the transition amplitude Eq. ( |3.7D maps out the zero-mode shape |'i/'o(3;)^| 
and the expectation value increases proportional to \fT for large T, 

(e-/'^-^±) ^/4^|Vo(x)p. (3.8) 

The case with a zero mode is numerically challenging, as is discussed in detail below. 

In the following sections, we present some exact expressions for (e^-/i' "^"^^i) for cer- 
tain inhomogeneous condensates ct(x), and we then compare these exact expressions to 
the numerical worldline results. These inhomogeneous condensates are taken from known 
inhomogeneous solutions to the gap equation of the GN model [50, 28, ^|. This serves as a 



precise test and control of our numerical scheme. These exact results are derived from the 
corresponding exact resolvents, derived in Appendix ^ using the Gel'fand-Dik'ii equation 
( ]B.4 , |B.5| ), by an inverse Laplace transform: 



1 ^oa 

i?±(x;-A) = (x|— r\x)= dTe-^^{x\e-^^^\x) . (3.9) 

+ A Jq 

For example, with a homogeneous condensate (t{x) = m, V± = m? , so the resolvent follows 
trivially from (p^). We find R±{ 
by an inverse Laplace transform. 

3.2 Single kink condensate 



trivially from ( p.4| ). We find R±{x\ — A) = ^-^=5=^, which leads to the familiar result ( p^ ) 



Consider the inhomogeneous kink-like condensate, studied by inverse scattering in |5C|, 

o-(x) = A tanh(^x). (3.10) 
The associated effective Schrodinger potentials (|3.2|) are 



[ 

V± = I . (3.11) 

[A^{l-2sedi^{Ax)) 

The kink condensate (t(x), and the corresponding V±{x), are plotted in Fig. ||. The poten- 
tial V-{x) has an exact zero mode, localized on the kink. 

The diagonal resolvents R±{x; — A) are derived in Appendix and an inverse Laplace 
transform gives the corresponding exact heat kernels: 



K±ix;T) 



/47rT 



MttT 



+ ^ Erf (aVt) sech2(Ax 



(3.12) 



where the error function is Erf(z) = e"* dt. According to Eqs. (^), (^, the 

zero mode of the V- potential is reflected in the large T behavior of K^{x;T): K^{x;T — > 
00) ^sech^{Ax) = |'(/'o(a;)p. 

Capturing the zero-mode physics is indeed a precarious if not pathological problem for 
the standard worldline algorithm discussed above. Close to the origin, the potential V- is 
negative, see Fig. ||, and consequently the exponent in e~ d-rV-ixir)) -g positive for certain 
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Figure 1: Plots of the kink condensate a{x) in Eq. (3.10) [solid, dark blue, line], and the potentials 
V±{x) in Eq. (3.11). The potential V- [dashed, red, line] is localized on the kink, while V+ [dotted, 
light blue, line] is constant. These plots are for A= 1. 



loops. Comparatively small loops with many points close to the origin can exponentially 
dominate the expectation value for large T, inducing an overlap problem for the Monte 
Carlo estimate. A similar overlap problem has been studied in [ [47[ . 

According to Eq. ( |3.8| ), the zero mode leads to a VT increase of the exact worldline 
expectation value. It is immediately clear that a finite and discrete loop ensemble cannot 
capture this property for T — > oo: since the distance between two neighboring points of a 
discretized loop scales with VT, any loop ensemble will eventually no longer resolve the 
negative peak of V- for large T, losing the information about the zero mode. Thus, for 
any finite, discrete loop ensemble generated with respect to the free worldline action, the 
expectation value eventually decreases exponentially in the large T limit. 

In Fig. ||, the analytic result for the expectation value {e~ 'lo d.^^-i^i'^))'^ ^t the origin, 
i.e., the center of the kink, for ^ = 1 is compared to the worldline numerical result for two 
loop ensembles with different numbers of loops. For small propertimes, the results agree 
nicely. However, for large T values, the numerical results indeed decrease exponentially; 
the values obtained using a smaller loop ensemble decrease more rapidly than the result 
computed with more loops. In this large T region, the statistical errors are much smaller 
than the real error, which is typical for an overlap problem. Let us stress, however, that 
the worldline estimate is still reliable for a much larger range of propertime than the 
first- or second-order heat-kernel expansion, which is an asymptotic small-T expansion [see 
Appendix B.5] that fails for propertimes larger than 0{1). These heat-kernel expansion 
approximations are also plotted in Fig. and we see clearly that the worldline result is 
far superior at large propertime T. 

We conclude that the zero-mode-induced overlap problem inhibits a straightforward 
calculation of the heat kernel's large T behavior and can only be shifted to larger T values 
by using a larger loop ensemble. In Subsect. we demonstrate that the overlap problem 
can be solved by generating an adapted ensemble with a Hybrid Monte Carlo algorithm. 
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Figure 2: Left panel: the CP-loop expectation value at the centre of a kink potential with A = 1 
versus the propertime T. For large propertimes, the analytic result is proportional to Vt. This is 
rediscovered by worldline numerics only up to a certain propertime value, depending on the number 
of loops employed. Still, the worldline estimate is reliable over a much larger range of propertimes 
than the first- and second-order heat-kernel expansion. Right panel: zero-mode wave function of 
the kink with A — 1 reconstructed from the worldline computation of the heat kernel at large T 
using a CP loop ensemble. 

In the present kink case, the zero-mode physics can still be identified with the standard 
free ensemble, since it is clearly separated from the remaining spectrum. Already at finite 
but large T, the heat kernel is essentially determined by the zero-mode contribution, with 
all other modes being exponentially suppressed. The right panel of figure |2| shows the zero- 
mode shape extracted from the heat kernel at a range of propertime values T of 0(10) using 
50000 CP loops consisting of 100 points per loop. In this propertime range, the modulus 
of the zero-mode wavefunction iV'oP can be reconstructed from ( |3.7D by a fit. Comparison 
with the exact result reveals that a rather rough discretisation of the loops already yields 
quite accurate results. 

3.3 Kink-antikink pair condensate 

The single-kink condensate discussed in the previous subsection has an exact zero mode 
for V^{x), which complicates the worldline numerics at large T. In this section, we probe 
this phenomenon more precisely by considering a condensate for which the bound state is 
shifted away from zero energy. Consider the kink-antikink configuration: 

cr(x) = A{coth{b) + [tanh(^a;) - tanh(^x -h b)]} 

= A{ [coth(6) - tanh(6)] + tanh(6) tanh(^ x) tanh(^ x + b)} . (3.13) 

which has the form of a kink-antikink pair, with finite separation b/A. The corresponding 
potentials ( |3.2D are 
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Figure 3: Plots of the kink condensate a{x) in Eq. ( t3.13| ) [solid, dark blue, lines], and the potentials 



V±{x) in Eq. ( 3.14 ). The first plot is for b = 0.5, while the second plot is for b = 10. Both plots 
have A = 1. Note that, in each case, the potential V- [dashed, red, line] is localized on the kink, 
while V+ [dotted, light blue, line] is localized on the anti-kink. 

The kink-antikink condensate (t{x), and the corresponding V±{x), are plotted in Figure |3|. 
The potentials in Eq. ( |3.14 ) are special in the sense that they are reflectionless, and each 



has a single bound state, located at ^^/sinh^(6), and a continuum starting at coth^(6). 
Note that in the hmit 6 ^ oo, the antikink disappears to minus infinity and we are left 



with the single kink configuration (t{x) = Atanh{Ax), as in Eq. ( 3.10 ). To see this 
pictorially, compare the second plot in Figure in the vicinity of x = 0, with Figure 
|l[ Correspondingly, the bound state at A^/sinh^(6) becomes a zero mode in this infinite 
separation limit. 

The resolvents R±{x; —X) for the potentials V± are derived in Appendix and the 
heat kernels follow from the inverse Laplace transform {\i.9\): 



. , I ........ (.^) l^^^l^l ^ ,3..) 

For finite b, the bound state at ^^/sinh^(6) is reflected in the exponential factor in the 
second term in (|3.15| ). Thus, the large-T behaviour isolates the lowest energy mode, the 
bound state at energy -Eb.s.^ 

i^±(x;r) ^ e-^'^-^|V'±'-(2:)|2 , T ^ oo. (3.16) 



In the limit b — > oo, this bound state becomes a zero mode, and the heat kernels in ( 3.15| ) 
reduce to those of the single- kink ( |3.12| ). In the limit b —>■ 0, A ^ 0, with A coth6 = m 
fixed, a{x) — > m, and K±{x;T) reduce to the heat kernels ( |3.6| ) appropriate for a homoge- 
neous condensate. For a comparison with worldline Monte Carlo results, we only consider 
the challenging bound-state contribution contained, e.g., in K-{x;T) [clearly K^{x;T) 
yields the analogous result for £+]. This corresponds to the second term in Eq. ( |3.15| ), 
yielding a bound-state contribution to the action density, i.e., the Lagrangian, 

£_,b...(x) = ^^sech^Ax) r ^ e-^^^/^-i^'W Erf(^^) - c.t. (3.17) 
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For illustration, we fix the parameter b by Acoihb = 1, so that the bound state energy 
is at -Bb.s. = 1 — A^. This choice minimizes the contribution of the continuum part of the 
spectrum to the action. Including explicit counter terms, satisfying Coleman- Weinberg 
renormalization conditions, we obtain (setting m = 1) 

/:_,b.s.(x) = ^ sech2(Ax) / ^ iy^e-^^-^"^'^ Erf(AVr) + 2A{e-^ - 1)) . (3.18) 

J 

The remaining free amplitude parameter A determines the depth of the V- potential: 
A = 1 corresponds to the single kink, with an exact zero mode, whereas ^ = yields a 
constant, positive potential, with no bound state; intermediate values interpolate between 
both extremes. If A is larger than 1/^/2, the potential becomes negative at the origin. But 
only for A = 1, there is a zero mode and the integrand in Eq. ( p. 18 ) is proportional to Vt 
in the limit T ^ oo. 

We have performed a worldline numerical computation for five potentials with A values 
between 1 / \/2 and 1, using the standard free ensemble in order to study the overlap problem 
quantitatively. The result corresponding to the propertime integral in Eq. ( 3.1^ ) is shown 



in Figure ^ This analysis shows that, in general, the negative minimum of the potential 
is not a problem. Only for the largest A value, A = 0.995, for which -Eb.s. = 0.00099975, 
does the algorithm have a severe convergence problem. In this case, the influence of the 
near-zero mode becomes noticeable: the exponential factor in the integrand of Eq. ( 3.1^ ) 



decreases only weakly, and the integral is dominated by large propertime values, at which 
the worldline numerical evaluation of the expectation value suffers from the overlap prob- 
lem, similar to that displayed in Fig. |2|. Whether or not a severe overlap problem occurs 
can be estimated from the near-zero mode properties. A near-zero mode creates a VT 
increase of the worldline expectation value at intermediate T values, being followed by an 
exponential decrease at very large T values. The crossover occurs at T values T ~ 1/E, 
where E corresponds to the near-zero mode's energy level. For instance, for the above 
kink-antikink case with A = 0.9995, the crossover occurs at T ~ (1 — A'^)~^ ~ 10^, imply- 
ing that a reliable estimate of the integrand up to values of T ~ 10^ would be needed in 
order to solve the overlap problem by brute force. For the case with A = 0.942809, the 
crossover occurs already at T~ {1 — A'^)~^ ~9 which is well resolvable by the standard 
algorithm with 1000 loops or more, cf. Fig. ^. 

3.4 Periodic array of kink-antikink pairs 

Next, consider a crystalline condensate formed from a periodic array of kink-antikink pairs: 

cn(6; i^)dn(6; v) 



a{x) 



sn o; u 



+ [Z{Ax]u) - Z{Ax + h;u)] 



A ( ""^ u sn(6; v) sn{A x; u) sn(^ x + b-u)] . (3.19) 

sn(o; v) 



Here Z{x]v) is the Jacobi zeta function, and sn, cn and dn the Jacobi elliptic functions. 



all with real elliptic parameter < < 1 |54, 55]. This type of inhomogeneous condensate 
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Figure 4: Convergence behaviour of the worldline-numcrical result for near-zero mode contribu- 
tions. The plot shows the numerical values corresponding to the propertime integral Eq. (3.18) 
versus the number of loops employed for different A values. The topmost values (green circles) are 
for the same A value as the values immediately below (solid grey squares), but with the number of 
loops multiplied by a factor of 1000. The dotted horizontal lines represent the analytic values. The 
error bars are the jackknife estimates and should provide reasonable values if the number of loops 
is larger than 1000. Evidently, this is not the case for the largest A value. This indicates an overlap 
problem. At about 2 750 000 loops {n — 2750), we observe a huge jump of the result, apparently 
caused by the occurrence of one single very small loop. 



is physically important in a number of contexts: it represents a bipolaron crystal in con- 
ducting polymers [56|, and it characterizes a crystalline phase in the massive GN model 
1 28]. In the limit ^ 1 we recover the kink-antikink configuration ( |3.13D , since in this 
limit: sn(5; z/) tanh(6), Z{h\v) — > tanh(6), cn(6; z^) — > sech(6), and dn(6; i/) sech(6). 
Note that we can restrict to < 6 < K(z^), because of the periodicity of the Jacobi elliptic 
functions. When b = K(z/), this condensate is important for soliton lattices in polymers 
1 56, for the Peierls model |5S], for periodic phases in superconductors [59|, and for the 



crystal phase of the massless GN model [28|. In this case, discussed in more detail in the 
next subsection, the neighbouring antikinks lie at the very edge of the period, and so in 
the infinite period limit \v ^ 1], the antikinks disappear to plus/minus infinity, and we are 

!T9|) contains all 



left with the single-kink configuration (3.1C). Thus, the configuration 
the other examples as special cases in a particular parametric limit. 

The potentials ( p.2| ) corresponding to the periodic condensate ( p. 19 ) are 



V± = A^ 



1 



sn2(6; v) 



l + v-2v 



cn^(^x -I- 6; v) 
cv?{Ax] v) 



(3.20) 
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Figure 5: Plots of the periodic condensate a{x) in Eq. ( 3.19 ) and of the periodic potentials V±{x 



in Eq. ( |3.20 ). The first plot is for elliptic parameter v — 0.5, for which K(0.5) = 1.85, while the 



second plot is for v = 0.999, for which K(0.999) = 4.84. In each plot, A = 1 and 6 = 0.5. The sohd 
[dark blue] curves are the periodic condensates. Note that V-{x) [dashed, red, lines] are localized 
on the kink, while V+{x) [dotted, light blue, lines] are localized on the antikink. Notice that as 
1/ ^ 1, the period diverges, and we recover the kink-antikink configuration: compare the second 
plot with the first plot in Figure ||. 

The kink-antikink crystal condensate cr(x), and the corresponding V±{x), are plotted in 
Figure ^. 

These potentials are periodic functions, with period 2K(z^)/^, and are special in the 
sense that they possess only a single bound band - they are the simplest example of the 
"finite-gap" potentials |Q, which are the periodic analogue of reflectionless potentials. 
V±{x) are the same function, displaced in x, and so they have precisely the same band 
spectrum. In the language of supersymmetric quantum mechanics, they are self-isospectral 
[61]. The band spectrum is plotted in Figure ^. The band edges lie at cs^ {b; v) and 
ds^(6; u), and with a continuum starting at ns^(6; u). As — > 1, the period diverges, 
and the bound band shrinks to a single bound state at ^^/sinh^(6), which is just the bound 
state of the reflectionless system associated with the single kink-antikink pair. 

The resolvents for the potentials in ( |3.20| ) are derived in Appendix and the corre- 
sponding propagators are obtained by an inverse Laplace transform: 



K^(.-,T)=c.(T)+mr'^t,T''\'\ ■ p-^i) 

cn[Ax;iy) 



Note that the coefficients a{T) and P{T) are the same for the two propagators - the only 
difference between and K- is the shift of the x dependence in the cn^ functions. The 
coefficients are most usefully expressed in the convolution form: 

/3(r) = ^ ^-A-^^H^:.)T /■ ^ / j„ ( (3.22) 



y^7r(r - U 



2 



a(T) = - { dsHb; v) m + ^) • (3-23) 
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Figure 6: The spectrum of the Schrodinger operators H± with (j{x) given by the periodic conden- 



sate Eq. (3.19), as a function of the eUiptic parameter v. The spectrum has a single bound band, 
with band edges A^cs^{b; v) and A^ds^(6; v), and a continuum band starting at ^^ns^(6; v). When 
I' — > 1, the bound band contracts to a single bound level at / sinh^ 6, with the continuum edge 
at coth^ b. In the other limit, as v ^ the spectrum is pure continuum, with lower edge at 
cot^ h. This plot is for = 1. These features can be clearly seen in the associated resolvents in 
the text. 

When = 1 we recover the propagators for the kink-antikink pair configuration in ( |3.15| ), 
while when = we find propagators for a homogeneous condensate of mass ^cot(6). 

An interesting special value of v is the lemniscate case, = ^, for which /3(T) takes 
the simple Bessel function form 



For a comparison with worldline numerics, we choose typical potential parameters as they 
occur in the phase diagram of the GN model at finite densities and temperature: 6=1, 
u = 0.9, and A = m/coth(l). Since this typical form of the potential is free of (near- 
)zero modes, it represents a direct test of the quantitative reliability of worldline Monte 
Carlo for fermionic fluctuations in generic backgrounds. Figure |^ shows the expectation 
value (exp(— j V-)) versus the propertime in a minimum of the potential V-. Due to the 
absence of a zero mode, the values decrease exponentially already for small propertimes, 
in contrast to the corresponding single kink result in Figure ^, and we see in Figure ^ 
that the numerical and analytic results agree very well. Even a small loop ensemble, of 
just 2000 loops, yields tiny errors. By contrast, the first- and second-order heat-kernel 
expansion, which is a standard analytic expansion technique, fails badly for T > 1. For 
further quantitative comparison, we integrate the resulting Lagrangian £_ over one period 
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Figure 7: The worldline expect at ion- value versus the propertime in a minimum of the potential 
V- for b — 1, V — 0.9, and A — m/coth(l). The statistical error is too tiny to be plotted. Due 
to the absence of a zero mode, the values decrease exponentially very early, in contrast to the 
corresponding values for a single kink, as shown in Fig. The figure above shows that the analytic 
heat-kernel expansion [the dashed and dotted lines] is not reliable for T > 1. 



along the x axis: 

r-2K(0.9)/A 



dx C- = — — 



(3.25) 



6.81704 (analytic) 
6.818 lb 0.001 (worldline numerics) 

where the upper value in braces is the analytic result, the lower one a worldline- numerical 
estimate with 20 000 CM loops of 2000 ppl each. (The result obtained for the integrated 
£+ is identical.) In units of m, the spatial average of C- then is 

2K{0.9)/k 



dx C- 



(3.26) 



2K(0.9) JO 
m2 J 0.766857 (analytic) 
471- \ 0.7669 ± 0.0001 (worldline numerics) 

Note that this is well above the corresponding value of the constant-field solution, — l/(47r), 
indicating that the constant field is the preferred vacuum solution. We conclude that a 
precision on the permille level or below is straightforwardly achievable with worldline Monte 
Carlo algorithms for fermionic systems in absence of (near-)zero modes. 

With the same loop ensemble we study the slightly more challenging configuration 
b = 2 , I' = 0.9 and m = j4coth(2), which has deeper minima, and obtain 

2K(0.9) ^ 

47r 



dx L- 



(analytic) 
(worldline numerics) 



0.499945 

'0 1^0.49 ±0.03 

with the spatial average (in units of m) 

TT? J 0.0901099 (analytic) 
4vr 1 0.088 ±0.005 (worldline numerics) 



(3.27) 



(3.28) 
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The statistical error is larger, as one might have expected, but can be improved using a 
larger worldline ensemble. The agreement with the analytic result is still very satisfactory. 



3.5 Kink crystal condensate 

An important special case corresponds to choosing the kink-antkink separation to coincide 
with half the period of the crystal, in which case the antikinks neighboring a given kink lie 
at the edges of the period. This corresponds to choosing b = K(i/). Then the condensate 
in Eq. ( 3.19| ) simplifies to 

a(x)=A. ^"^^;-;^^"^^"'"^ (3.29) 
dn(^ x; v) 

The associated Schrodinger potentials are 



(1 - u)sd^{Ax]u) 
cv?{Ax; v) 

The condensate and associated potentials are plotted in Figure 



V. = A^\.-2.r^ — — --'H . (3.30) 




Figure 8: Plots of the periodic condensate ij{x) [solid, dark blue, lines] in Eq. (3.29), and the 
periodic potentials V±{x) in Eq. ( [3.30| ), with elliptic parameter v = 0.5 [first plot] and v = 0.999 
[second plot]. Both plots are for ^ = 1. Note that the antikinks lie at the edge of the period, 
so that V^{x) [dashed, red, line] is localized on the kinks at the center of a period, while Vj^(x) 
[dotted, light blue, line] is localized on the antikinks at the edge of a given period. 

This system has a band spectrum with band edges at 0, A^{\ — v), and A^, as shown 
in Figure |9|. Note that the bottom of the lowest band is at 0, for all and that there 
is just one bound band. The corresponding heat kernels K±{x]T) can be obtained from 
those in the previous section by the simple substitution b K(z^). 

3.6 Worldline ensembles from Hybrid Monte Carlo 

Since the standard worldline ensembles generated with respect to the free worldline action 
can lead to an overlap problem for the special case of a (near-) zero mode in the fermion 
spectrum, let us explore a different route: we generate the ensemble with respect to the full 
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Figure 9: The spectrum of the Schrodinger operators H± with (j{x) given by the periodic conden- 



sate Eq. (3.29), as a function of the elliptic parameter v. The spectrum has a single bound band, 
with band edges and A^{1 — i'), and a continuum band starting at A^. When v ^ 1, the bound 
band contracts to a single bound level at 0. In the other limit, as — )■ 0, the spectrum is pure 
continuum, with lower edge at 0. 

action using a Hybrid Monte Carlo (HMC) algorithm. In principle, worldline generation 
with the full action would also be possible with a local update algorithm such as the 
Metropolis update. However, due to the fact that the worldlines are one-dimensional, local 
updates suffer from a long auto-correlation time, as observed in [^]. This is avoided by 
the global HMC update. 

The quantity which will be central for the HMC simulation is given by 

K±{x,T) = (x|e-^(-^?+^±)|x) . (3.32) 

This quantity is given by means of an expectation value with respect to the full loop 
ensemble: 

W{x,T) = /-dl + V±) , (3.33) 



xF 

where the loops are generated with the probability 

PfN = 5 [xcM - dTx{T)^ exp ^ dri;2(r) - ^ dTV±{x{T)) . (3.34) 

The strategy is to generate statistically important loop configurations {x(r)} using HMC 
methods and to estimate W{T,x) in ( 3.33| ). The advantage of this approach is that the 



loop ensemble already incorporates features of the potentials V±, and that far fewer config- 
urations are required for a reliable estimate of expectation values than for free loop clouds. 



-19- 



The disadvantage is that more computational resources are necessary to generate the HMC 
loop ensembles. Ultimately, it will depend on the particular potential whether the addi- 
tional numerical effort is justified. We will argue that this is the case for a Hamiltonian 
H± = {—df + V±) which possesses zero modes. 

Once a numerical estimate for W{x, T) is obtained, the heat kernel K± can be re- 
constructed by integration of ( 3.31| ) with respect to the propertime. Note that the kernel 
satisfies the boundary condition K±{x,0) = 1. We therefore find: 

K±ix,T) = expj- ^ dTW{x,T) | . (3.35) 

Let us briefly comment on the implementation of the discretized loops. As before, the 
propertime interval r = [0, T] is represented using points: 

T 

T = ndr , dr = — , n = l...N , N odd . (3.36) 

For a simplification of the notation, we will assume that is odd. We base the HMC 
algorithm on a discrete Fourier representation of the worldlines, which can easily account 
for the worldline periodicity and leads to uncoupled degrees of freedom for V± = 0. Hence, 
considering the Fourier coefficients as degrees of freedom implies that loop ensembles can 
be generated with little autocorrelation as long as the influence of the potential is not too 
significant.^ For CM loops, we have the representation: 



JV-1 
2 



x{t) =xcm + -^ ^ Cfcexpfi^fcrj , c_fc = 4, cq = 0. (3.37) 
In the case of CP loops, for which xcp is the common point, we analogously find 



JV-1 
2 



x{t} = Xcp + ^ J2 '^f' 



(3.38) 



With these representations, the worldline kinetic term is GauBian in the weights c^: 

4 



n=l k 

The task is now to generate statistical ensembles of the degrees of freedom 

k - -1 1 ^ 



which are distributed according to, e.g., the CM probability distribution ( 3.34 ). Note that 
the (5-function constraint has been exactly incorporated leaving us with 2 ^^^^^ = N — \ 



^An HMC algorithm could equally be based on the v loop algorithm since it features the same 
necessary properties. By contrast, an HMC generalization of the d loop algorithm j5l| which is based on 
an iterative loop construction seems not to be straightforward. 
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degrees of freedom, parameterized by the coefficients c^. For the HMC approach, it is highly 
advisable to introduce real valued degrees of freedom. We therefore define the coefficients 
dfc G R by 

N -I 

dk = Cfc + c_fc , d^k = i {ck - c^k) , k = l... — - — . (3.39) 

With this definition the loop coordinates, for instance, for the CM loops ( |3.37| ) are mani- 
festly real valued: 

N-l 

x(r) = xcM + ^ 1^ cos^^/cr^ + sin A; | . (3.40) 

The next step is to introduce the canonical momentum vr^ G M for each degree of freedom 
dk- The corresponding HMC Hamiltonian is defined by 



k 

V(4) = 7 r dTx\T) + r dTV±{x{T)) . (3.42) 
4 Jo Jo 



The degrees of freedom experience an evolution with the so-called HMC time u. This time 
evolution is determined by the Hamilton equations of motion: 

d , d dVjdk) 

— dfc = TTfc , -j-^k = — • (3-43) 

au du odk 

The important property of any solution dk{u), vrfc(u) is that the HMC energy is conserved: 

^ n(7rkiu),dk{uj^ = . 

The initial momenta vrfc(n = 0) are selected at random from a Gaussian distribution. The 
initial coordinates Ck{u = 0) are given by their actual ensemble values. With these initial 
conditions, the values at u = Uf are approximately obtained with the help of a Leapfrog 
integration of the equations of motion ( p. 43 ). Because the equations of motion are not 
solved exactly, the HMC energy 7i is not exactly conserved. The idea central to the HMC 
approach is to accept the values dk{uf) as new members of the Markov chain with the 
probability 



p = min(^l, exp{-AF}j, /\H = W ^^^(nj), 4(n/)J - ?^ ^71^(0), 4(0)J . (3.44) 

Since the equations of motion are satisfied to a good extent due to the Leapfrog integration, 
the value Ck{uf) can be largely different from Cfc(O) (implying good ergodicity) while the 
acceptance rate is still high. 

In order to test the HMC approach, we revisit the kink background field ( 3.10| ) (with 
A = 1 in the following). The HMC final time is given by = N^r&du where A'tra is 
the number of steps of the Leapfrog integration. The step size du is chosen to obtain an 
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Figure 10: Left panel: the HMC energy as a function of step size du of the Leapfrog integration. 
Right panel: the function W{x = 0, T) for the kink configuration obtained with HMC loops. The 



large-T behavior is dominated by the ground-state energy, of. Eq. (3.45). Sufficiently many points 
per loop N are required in order to observe the approach to the zero-energy mode of the kink. 



acceptance rate between 20% and 80%. Figure (left panel) shows the HMC energy as 
function of the step size du for several values for the proper time. We find that a value 
du PS 0.3 provides a good compromise between ergodicity and acceptance rate. Using the 
HMC ensembles for CM loops, the function W{x, T) is estimated and the heat kernel is 
reconstructed with the help of the proper time integration (see Eq. ( f3.35D ). Note that 
because of ( 3.31| ), we find for large values of T 

|2 



W{x,T) ^ Eo + (Ei-Eo) 



K^|l) 



,^-(E,-Eo)T large), 



(3.45) 



Kx|o)p 

where £"0, Ei are the energies of ground and first excited states, and |0) and |1) are the 
corresponding wave functions. In theories with a gap between the first excited and the 
ground state, the latter equation implies that W{x,T) rapidly approaches a constant, i.e., 
the ground state energy Eq, for large T. In the worst case scenario of a zero mode, W{x, T) 
approaches zero and a positive constant otherwise. Figure ^ (right panel) shows the re- 
sult for W{x = 0,T) for the kink configuration. Indeed, the function W{x,T) rapidly 
approaches a constant for large T. The overlap problem, discussed in the previous subsec- 
tions, has been solved by the use of HMC loop ensembles. Note, however, that, also in the 
present case, the number of points per loop N must be large enough: for coarse loops, the 
numerical estimate approaches a negative constant thus erroneously suggesting a negative 
ground state energy and an instability of the operator. A sufficiently fine representation 
or an improved integration over loop coordinates might resolve this problem. 

~ We 



Finally, we present our final result for the reconstructed heat kernel in Fig. 11 



obtain very good results with HMC using only a modest number of points per loop. 
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Figure 11: HMC results for the heat kernel for the kink configuration in comparison with the 
result obtained from free loop ensembles (all used ensembles are CM loops). Significantly fewer 
HMC loops are required for a reliable estimate of the large-propertime behavior. 

To conclude this subsection, the overlap problem of the free loop cloud method is 
indeed solved by adopting the HMC approach. Note, however, that a precise estimate of 
the effective action in the case of a (near-)zero mode is still a numerical challenge. 

4. Worldline approach at finite temperature and chemical potential 

Even though introducing finite temperature and chemical potential appears straightforward 
in the worldline approach, special attention has to be paid to combining this straightforward 
formalism with Monte Carlo algorithms. The reason is that an analysis of field theories at 
finite temperatures and densities often requires knowledge about the analytic structure of 
propagators etc. in the complex energy-momentum plane, or similar conjugate variables 
such as the propertime. By contrast, Monte Carlo algorithms require a positive real action 
for importance sampling, and this restricts the computable information to lie on the real 
axis of given variables. 

Since this conflict of interests between thermal field theory at finite density and Monte 
Carlo algorithms is general and not only restricted to the Gross-Neveu model, we discuss 
this problem in more detail in the following. 

4.1 Temperature in the worldline representation 

Consider the worldline formulation of a general quantum-field theoretic amplitude. Finite 
temperature can easily be implemented with the aid of the Matsubara formalism: The 
Euclidean time, say along the Z^th direction, is compactified to the interval [0, /?] with 
periodic boundary conditions for worldline fluctuations for bosonic fields and antiperiodic 
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boundary conditions for fermionic degrees of freedom. Here f3 = 1/T is the inverse tem- 
perature. (Note that we reserve the symbol T for propertime, whereas the temperature is 
denoted by T.) As a consequence, the worldhnes can also wind around the time dimension. 
It is convenient to write a given loop x(r) with winding number n as sum of a loop with 
no winding, x(r), and a translation in time running from zero to n/3 with constant speed, 

T 

x^,{t) = x^{t) + n(3-5^D- (4.1) 

The path integral over the different winding number sectors labeled by n factorizes for 
static configurations, yielding 



L 



Vx e- ^0 dr^- . . . = V (-l)"e-'^ / Vxe-^od^'r .... (4.2) 

x(0)=x(T) Jx{0)=HT) 

where the factor (—1)" implements antiperiodic boundary conditions for the present fermi- 
on fluctuations; it would be absent for bosons here and in the following. A finite-tempera- 
ture generalization of a vacuum propertime expression is straightforwardly obtained from 
a simple insertion into the propertime integral. 



(4.3) 



where the ellipsis represent, for instance, a heat-kernel expression or a worldline expectation 
value etc. The winding number sum over n is directly related to the Matsubara sum by a 
Poisson transformation. Quantum and thermal fluctuations are separated by the different 
winding number sectors: the n = term is identical to the vacuum expression, while the 
|n| > 1 sectors are purely thermal contributions which vanish in the vacuum limit /3 — > oo. 

Let us stress that the present thermal formalism not only holds in general but is also 
completely robust when combined with worldline Monte Carlo techniques, as has been 




used in the context of the Casimir effect at finite temperature [32|.. The winding num- 
ber sum converges rapidly in the relevant low- and intermediate temperature range; at 
high temperatures where increasingly more winding sectors would have to be resummed, a 
Poisson resummation to Matsubara sums can be performed, such that rapid convergence 
is guaranteed for all temperatures. In particular, analytical or numerical knowledge of 
the heat kernel on the real T axis is sufficient for a reliable estimate of thermodynamic 
amplitudes. For example, it is straightforward to derive the critical temperature for sym- 
metry restoration in the Gross-Neveu model within the present formalism, yielding the 
well-known result 

= Tre-'^m"^ PS 1.76387m- \ Tc = 1/^c ~ 0.566932m (4.4) 



in units of m |63|. For temperatures larger than the critical temperature 7^, the potential 



a vanishes and the spontaneously broken discrete chiral symmetry is restored. 
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4.2 Chemical potential in the worldline representation 

Introducing a finite chemical potential via the grand-canonical partition function leads us 
to the microscopic action of the Gross-Neveu model at finite density, 

S = Jd^'x (^-^ + iP + , (4.5) 

in D dimensions. The resulting contribution to the effective action from fermionic fluctu- 
ations reads analogously to Eq. ( ^ ) 

Ar[a] = -^Trln(-9^ - {Od + l^f + - i^a). (4.6) 

Since the chemical potential occurs in the form of the Dth component of an imaginary 
Euclidean gauge potential in the effective action, it is tempting to jump directly to the 
corresponding worldline representation in the presence of gauge fields, which are coupled 
to the worldline in the form of a Wegner- Wilson loop. In other words, the free worldline 
action receives a further contribution of the form 

g-i/o^'^™^ Q-\!o drx'^-t^Io '^■ri'D (_ i)"e~ e~ 3 /o^ (4.7) 

where in the last step we have used the decomposition of worldlines into different winding 
number sectors n on the finite-temperature cylinder. Since also negative winding numbers 
are allowed, it is already apparent that the winding number sum may suffer from conver- 
gence problems at large This issue will be discussed in detail in the next subsection; 
let us simply state here that convergence problems are absent for (3^ < vr which we will 
assume in the following. 

Decomposing the winding number sum into the vacuum sector n = 0, and the temper- 
ature and density corrections \n\ > 1, implying a decomposition of the effective action (or 
grand potential) into a vacuum part F, and a thermal part AT^^^: 

r^,^ = r + AT^,^, (4.8) 

with r given by Eq. (|2.18| ), we end up with the worldline representation for the thermal 
contribution at finite density 

(4.9) 

For Pfj, < vr, this representation is reasonably converging and can numerically be imple- 
mented rather straightforwardly. Naively interchanging summation and integration yields 
rapidly converging propertime integrals for any given n; however, the n sum, which is 
reminiscent of a fugacity expansion, then diverges for > vr. Of course, Eq. ( [4.9D can 
immediately be used for a finite-density expansion in to any order. It is also straightfor- 
wardly applicable to the case of imaginary chemical potential where all convergence issues 
are absent. 
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Figure 12: Phase diagram for a constant potential (cf. | |63[ ) computed by numerically minimising 
the full action with respect to a. The colour corresponds to the resulting a value, the phase of 
restored chiral symmetry is in black. The tricritical point is marked by a white circle, located 
between phase transitions of second order (the upper boundary of the blue region) and transitions 
of first order (discrete jump of the colour/cr value). The line of second order transitions extends to 
the /I = axis, where it is marked by an arrow at the value % ~ 0.567 Eq. (4.4). The blank region 
below the T = ji/ir line is inaccessible to the given expression for the action as discussed in the 
text. 



As an application of Eq. ( [4. 91 ) , let us briefly revisit the phase diagram for a homogeneous 
a condensate. The worldhne expectation value for this is trivial, since exp(— JV±) = 
exp(— Tcr^), for any single worldline. This case serves as a simple but instructive test 
for the (non- worldline) numerics of our algorithms, and for the representation Eq. ([4.S|). 
In Fig. |l2|, the phase diagram obtained by numericahy minimising the full action with 
respect to a for a given chemical potential and temperature is shown. In the accessible 



domain, that is for < vr, the diagram agrees with the result of [33|, as it should. At 
/i = 0, we recover the critical temperature in Eq. ( |4.4| ). Even close to the line = vr, i.e. 
temperature T = /x/vr, below which the propertime integral has convergence problems, we 
obtain a stable result. 

4.3 High densities and low temperatures 

Let us further investigate the instability for /3/i > vr. For this purpose, we first address the 
homogeneous background (T=constant, <I>[cr] = 1, and extend the findings to the general 
case at the end of this subsection. 

The (regularized) effective action is given by 
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In order to perform the sum over the windings n, we introduce a Hubbard-Stratonovich 
transformation by 



i-ire 



4T e 



-fi/Bn 



UttT 



/52 



dv e 



(4.11) 



Inserting (4.11) in ( [4.10 ) yields for the T integral 
1 f .n N,d^ 1"°° dT ' 



Nidy dT MvrT 



(4.12) 

The IR problem for large values T has become apparent: the T integration converges 
provided 



Re 



> 0. 



(4.13) 



For imaginary chemical potential this condition is always satisfied. But for real chemical 
potential we find the condition 

< \/^2 + p2^2_ (414) 

Allowing for the possibility of a vanishing condensate, or in the more general case for the 
existence of a zero mode in the spectrum, we restrict to chemical potentials satisfying the 
bound 

PfKTT. (4.15) 



The basic idea of the present subsection is to consider the case ( p^ ), then to der ive 
an analytic expression for the effective action, and to perform an analytic continuation to 
values ^ > 7r//3. For a homogeneous background, we find that this approach recovers the 
well-known result familiar from solid state physics. 

Hence, let us confine ourselves to values /i satisfying ( [4.15| ), and perform the T inte- 
gration first. For this purpose, we introduce the "density of states" po by 

^ - 2(47r)(^-i)/2 



2^(D-l)/2 



dEEpo{E) e 



Po{E) 



(27r)^-i 

Rearranging ( 4.12 ) with the help of ( [4.16 - 4.17 ) yields 

/roo 1 ^ poo jrp 

d'^X / dEEpo{E)- / — 



(4.16) 
(4.17) 



(4.18) 



The n sum, the v integral, and the proper time integration can be performed in closed 
form. The details of this calculation are presented in appendix]^. The final result is: 



-Am 



roo,o 

Nfd-y 



(4.19) 



j d^'x dEEpo{E) [in (l + e-^(^-^)) + In (l + e-^(^+'^)) 
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which we recognize as the standard expression [|6^ . Although this result was derived under 
the assumption ( 4.15| ), the closed form ( [4.19| ) is perfectly valid for > vr. Hence, we 



perform an analytic continuation by taking ( [4.19| ) as the primary expression and consider 
the arbitrary values for in the following. 

Let us generalize the above findings to the more general case of a spatially dependent 
background field: 



(4.20) 

The key observation is that we can relate the general case ( |4.20| ) to the free case by defining 
the density of states p{E) by 

j;^{e~^o'^'''m) = 2{A.r~'y' j^dEEp{E)e-^''\ (4.21) 

The final answer for the effective action F^^^ is given by (f4.19| ) where we replace /?o by p 
from ( [4.21 ). Note that the left-hand side of (4.21) is only known numerically for real values 



T. The quite challenging task is then to perform the inverse Laplace transform to obtain 
p{E). We however point out that, depending on the parameters, the full range of p{E) is 
not needed. Let us consider the particle density np^^ for the moment {V = j d^x), 



1 5r 

np,^ = -- = Nfd^ dE E p{E) 



1 



(4.22) 



For small temperatures, i.e., for (3p, ^ 1, we obtain the result ( valid if /x lies within a 
band) : 

np,t. = Nid^l^iyEEp{E) + ^±[Ep{E)]\^^^ + 0{l/(3^)Y (4.23) 

In this case, the density of states need only be known for E ^ p,. 
4.4 Crystalline phase on the worldline 



The crystalline phase of the GN model extends up to the tricritical point plotted in Fig. (12) 
p8[| . This point and therewith some part of the crystalline phase is above the line defined by 
Pp = TT, where Pp > vr. Consequently, already the straight forward evaluation of Eq. ( ^Ig| ) 
allows for studying translational-symmetry breaking. 

We consider several points on a line slightly above the Pp = vr line, given by 1.1/3/i = vr, 
which intersects the crystalline phase between p = 0.63m and p = 0.67m. For each point. 



we determine the potential a of the ground state analytically as described in [30|. Using 



this potential, we evaluate Eq. ( |2.18 ) and Eq. ( |4.9| ) with worldline numerics and compare 



the resulting value of the grand potential per volume with the corresponding value for the 
constant field solution. 

Figure (|^) shows the difference of these two values as well as the analytic result ac- 
cording to [50]. Both results agree nicely. In particular, the worldline numerical result 
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Figure 13: Difference between the grand potential per volume for the crystal and the homogeneous 
solution for T = I.I/x/tt. The solid line is the analytic result from js^, the crosses with error bars 
are worldline numerical values. The line T = I.I/x/tt crosses the crystal phase between /i = 0.63to 
and 0.67m, approximately. Consequently, we obtain negative values in this range. The numerical 
values agree with the analytic result and are significantly below zero. 

confirms the thermodynamic preference of the crystal solution over the homogeneous po- 
tential in the range 0.63m < fi < 0.67m which is reflected in negative values of the plotted 
difference. 



5. Conclusions 

We have demonstrated that the Monte Carlo approach to the worldline formulation of quan- 
tum field theory may be consistently implemented to study interacting fermionic models 
such as the Gross-Neveu model. We have explained the required formalism and tested 
numerically the worldline results in comparison to a number of analytic results for inho- 
mogeneous condensates, pointing out the advantages and disadvantages of the numerical 



approach. In particular, in Section 4.4 we have used this approach to confirm the existence 
of the inhomogeneous crystalline phase of the Gross-Neveu model. While this is a highly 
non-trivial test, the true potential of the worldline Monte Carlo approach lies in the pos- 
sibility of scaling up to higher dimensions without dramatic numerical problems, and the 
possibility of realizing chiral symmetries directly. These issues are under investigation. 
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A. Fermionic models and continuous chiral symmetry 

In this work, we have mainly focussed on the the D = 1 + 1 dimensional Gross-Neveu (GN) 
model defined by the Euclidean fermionic action in (|2.1| ). The GN model has a discrete 
chiral symmetry under ip — > 75V', rather than a continuous chiral symmetry. In order to 
discuss the fate of chiral symmetry in our numerical worldline approach, it is instructive 
to study the Thirring model, defined by the action 



5f = y + ^(V'TmV')' j , (A.l) 

In addition to global phase and axial rotations, the Thirring model is invariant under 
SU(A^f)RX SU(A^f)L chiral transformations. By means of a Hubbard-Stratonovich trans- 
formation, the model can be rewritten with the aid of a vector boson 



'S'fb = / (fx 



2g^ > 



(A.2) 



where is invariant under chiral transformations. In the large- Af limit, the effective 
action upon integrating out the fermion fiuctuations yields 

nVf^] = I d^'x^y^ - Af Trln(-^ - if), A^f ^ 00. (A.3) 

As in the Gross-Neveu case, the resulting Dirac operator P = — f has a 75 hermiticity, 
2?t = ry^D^^^ implying that F can be written as 



T[V^] = / d'^x-^V^ - -iTVln -D[Vf + -a^,V^, , Af ^ 00, (A.4) 



^^-^V^ - ^TVln (-D[Vf + -u,,.., 

where D\y\ = + iV^, a^u = and the vector field strength V^,^ = d^Vy - duV^i- 

By virtue of the D = 1 -|- 1 dimensional Dirac algebra, the spin-field coupling (Pauli term) 
can also be expressed as (Tf^^V^i, = ^^e^yV^y = "f^V. Here, we prefer to work with the 
representation given in Eq. ( [A.4| ) which generalizes straightforwardly to higher dimensions^. 
Note that the fermionic contribution exhibits an accidental local gauge symmetry, which 
is, however, not respected by the classical action ~ V^. 

Incidentally, the 4-fermion vector interaction of the Thirring model is the only SU(Af)R 
X SU(Af )l chirally invariant and parity-even interaction which contributes to leading-order 

^Odd dimensions require reducible representations of the Dirac algebra in order to maintain 75 hermitic- 
ity of the Dirac operator. 
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at large- A'^f in D = 2. In higher dimensions, also independent axial-vector interactions exist 
which, in Z) = 2, are isomorphic to the vector interaction. 

The worldline representation for the Thirring model is well known from the analogous 
QED case. Concentrating on the contributions from the fermion fluctuations, we obtain 

AT[V] := -^Trln (^-D[Vf + ^a^uF^, 

oo 

= f(4;^'-/rW / ^--Uo^-^V)e-^f-.n.<,[y], (A.5) 

1/A2 x(0)=x(T) 

where the spin factor this time reads 

^V] = ^tr^ Vexp (^-\ dTG^^V^^ . (A.6) 

B. Exact resolvents and propagators for (1 + l)-dim. Gross-Neveu model 

In this Appendix, we present the derivation of the exact expressions for (e-io 'i^v±) dis- 
cussed in Section ^ where they were also compared to our numerical worldline results. We 
compute, for various forms of the static condensate function ct(x), the worldline propagators 

\/47rT 

where the Schrodinger-like Hamiltonian operators H± are 

H^^-dl + V^{x) = -dl + a\x)±a'{x) . (B.2) 

Our strategy is to consider the associated "diagonal resolvents" [i.e., coincident-point 
Green's functions], which are Laplace transforms of the propagators: 



R±{x;-\) = {x\— r\x)= dTe-^^{x\e-^^^\x) . (B.3) 

^± + A Jo 

For any static condensate function a{x), the resolvent must satisfy the Gel'fand-Dik'ii 
equation [S5| 

-2R R" + {R'f + AR^ (V + X) = 1 , (B.4) 



or its third-order linear form [obtained by differentiating ( |B.4 )] 

R'" - AR' (V + X)- 2RV' = (B.5) 

This identity follows directly from the fact that the Green's function in 1 dimension can 
be written in terms of the product of two independent solutions, and the product of two 
solutions necessarily satisfies such an equation. For details, see [35|. 

For certain special forms of static condensate cr{x), these Gel'fand-Dik'ii equations 
are analytically soluble. Then K±{x;T) is obtained as the inverse Laplace transform of 
R±{T; —A). Let us consider several examples: 
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B.l Homogeneous condensate 

With a homogeneous condensate cr(x) = m, we have V± = m^, so that solving (|B.4| )-(py 
leads to 

R^{x;-X) = -y2=== . (B.6) 

It is straightforward to take the inverse Laplace transform to arrive at the free propagator 
in (13^ 



B.2 Single kink condensate 

Consider the inhomogeneous kink-like condensate (t{x) = A tanh(^x), for which 

V± = < (B.7) 
[A^{l-2sedi^{Ax)). 



Then R±{x; — A) can be obtained from ( |B.4| )-( py5| ) by substituting the ansatz form 

R±{x;-X) = a±{\) + P±{X)sech'^{Ax) , (B.8) 

and solving algebraically for the [A-dependent] coefficients a± and f3± . A simple calculation 
leads to 



R±{x; —A) = < 



The corresponding propagators follow from the inverse Laplace transform, and are pre- 
sented in 

B.3 Kink-antikink pair condensate 

Consider the inhomogeneous condensate corresponding to a kink-antikink pair 

cj(x) = A (coth(6) + [tanh(^3;) - tanh(Ax b)]) 

= ^([coth(6) -tanh(6)] -htanh(6)tanh(^x)tanh(Ax-F6)) . (B.IO) 

The corresponding potentials are 



The resolvents can be derived by similar ansatze to ( [B.8[) , using sech^ {Ax + b) for V+ , and 
sech.'^{Ax) for V^. This leads straightforwardly to 

-A) = ^ ' ( 1 + H'lt.: \'^] I (B^12) 



The propagators again follow from the inverse Laplace transform, and are presented in 
( |3l5D - 
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B.4 Periodic array of kink-antikink pairs 

Consider the condensate formed by a periodic array of kink-antikink pairs 

cn(6; i^)dn(6; z^) 



a{x) = A 



sn o; V 



+ [Z{Ax]u) - Z{Ax + b]u)\ 



A ( — ^-^—r, — -r^—^ + V sn(6; v) sn(^ x; v) sn(A x + h\v^ 



sn(6; v) 

The potentials corresponding to the periodic condensate ( |B.13D are 

1 



V^ = A' 



sn^(6; v) 



cn^(^x + 6; v) 



(B.13) 



(B.14) 



The resolvents for V± can be obtained by substituting ansatze as in (|B.8| ), but with sech^ 
replaced by cn^, and with the argument shifted by h in the case of V+. Straightforward 
algebra yields 



R±(x--\) 



1 + 



A^v 



CY?{Ax + 6; v) 



(B.15) 



The propagators again follow from the inverse Laplace transform, and are presented in 
p.21 ). The kink crystal case, where the antikinks lie at the edges of the period, is obtained 
by the simple substitution h K(i/). 

B.5 Comparison w^ith heat kernel expansion 

It is instructive to compare the exact transition amplitudes, or heat kernels, not only with 
our numerical worldline loop results, but also with the approximate small T heat kernel 
expansions [l66| , 67, 38 1. In general, the transition amplitude/heat kernel has an asymptotic 
small T expansion 



^ oo 
n=0 



(B.16) 



where the functionals bn{x) are expressed in terms of the potential and its derivatives, and 
are given by: 



1 



-V 



6o(x) 
bi{x) 



b3{x) 



(B.17) 
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Defining 6„(x) 



2n + l 



(2n-l) 



jYTT r„(x), the r„ are given recursively by : 



^ n— 2 ^ 71—2 n— 1 n— 1 



fc=0 



fc=l 



fc=0 



k=l 



with T-Q = I and ri = —\V. 

For example, for the single-kink condensate, with V- = (l — 2sech^(74x)), the heat 
kernel expansion (B.16) is 



K_{x;T) 



i^l-A^T [1 - 2sech2(Ax)] + \a't' 



1 — - sech^(Ax) 



6 



6 

1 sech2(^x) 

5 



(B.19) 



in agreement with the small T expansion of the exact result in ( p. 121 ). This comparison 
is made in the first plot of Figure 0. Similarly, for the periodic kink-antikink condensate, 



with V-{x) as in (B.14), the heat kernel expansion (pTTq) is 



K_{x-T) 



1 



^/A■nT 



1- A^T 



U-1 + 



1 



sn2(6; u) 



2 cn^ (A x; v) 



+-A^T^ 



(l-z.)(^ + 3) ^ 2{v-l) ^ 1 



3 
4z^ 



sn2(6;z^) sn^(6;z^) 
v{v + l) \ cn^ {A X\ V) 



+ ... 



(B.20) 



in agreement with the small T expansion of the exact result in (3.21). This comparison is 
made in Figure ^. 



C. Effective action for the homogeneous case 



This section provides detailed information on the evaluation of the integrals and sums 
which appear in the effective action ( 4.1^ ) for a homogeneous background field. We firstly 
focus on the contributions from the terms with n 7^ 0: 



(C.l) 



Let us first perform the propertime integration. Using the identity 



2 / dxx e-^-' = le-^^' 



(C.2) 



the T integration can be be performed: 



(C.3) 
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The V integration requires to distinguish the cases x> \jl and x < //, but can 
straightforwardly done using the Cauchy integral theorem: 



otherwise be 



1{E) = dx e-^"(^+'') + sign(x - /x) e" 



/3n\x~iJ,\ 



(C.4) 



n=l 

The geometric series in n can easily be summed: 



m 



Using the identity 



the expression 



I{E) 



dx 



1 _|_ g-/3{a;+At) 



+ sign(x — fi) 



1 _|_ g-/3|x-/i| 



(C.5) 



-f3{fM-x) 



can be written as 



dx 



E 



-P(x+iJi) 



g-/3(x-/i) 
1 _)_ Q-l3(x-fj.) ' 



9{n-x) + 



-0{x-tM) 



Let us now consider the contribution with zero winding n = in Eq. ([4.18| ): 



(C.6) 



\ r°° dT f 



(C.7) 



It is convenient to decompose the latter expression as 



Io{E,fi) = AIoiE,fi)+IoiE,0), 



AIo{E,fi) = -l —e 



dT 



dv 



(C.8) 
(C.9) 



Note that the latter expression is UV finite. We have therefore removed the regulator by 
taking the limit A — > oo. A direct calculation of AlQ{E,fi) along the same lines outlined 
above yields 



AIo{E,i^) 

Thus, combining (|C.6|) and (C.7) gives 



dx 6{fi — x) 



(C.IO) 



I{E)+Io{E,fi)=Io{E,0)- I dx 

IE 



^-I3{x+fj.) g- 
+ 



■I3(x-^i) 



1 _)_ g-/3(x+/x) ]_ _|_ g-/3(x-/x) 



(C.ll) 



The final integration can easily be done leaving us with 

1 



I{E)+Io{E,^i)=h{E,Q)-- 



ln(l + e-^(^-'^) +lnfl + e-^(^+^) 



(C.12) 
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